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Abstract 

We perform a linear dynamical stability analysis of a general hydrodynamic model of 
chemotactic aggregation [Chavanis & Sire, Physica A, 384, 199 (2007)]. Specifically, we 
study the stability of an infinite and homogeneous distribution of cells against "chemo- 
tactic collapse" . We discuss the analogy between the chemotactic collapse of biological 
populations and the gravitational collapse (Jeans instability) of self-gravitating systems. 
Our hydrodynamic model involves a pressure force which can take into account several 
effects like anomalous diffusion or the fact that the organisms cannot interpenetrate. We 
also take into account the degradation of the chemical which leads to a shielding of the in- 
teraction like for a Yukawa potential. Finally, our hydrodynamic model involves a friction 
force which quantifies the importance of inertial effects. In the strong friction limit, we 
obtain a generalized Keller-Segel model similar to the generalized Smoluchowski-Poisson 
system describing self-gravitating Langevin particles. For small frictions, we obtain a 
hydrodynamic model of chemotaxis similar to the Euler-Poisson system describing a self- 
gravitating barotropic gas. We show that an infinite and homogeneous distribution of cells 
is unstable against chemotactic collapse when the "velocity of sound" in the medium is 
smaller than a critical value. We study in detail the linear development of the instability 
and determine the range of unstable wavelengths, the growth rate of the unstable modes 
and the damping rate, or the pulsation frequency, of the stable modes as a function of the 
friction parameter and shielding length. For specific equations of state, we express the 
stability criterion in terms of the density of cells. 

Key words: Nonlinear mean field Fokker-Planck equations, generalized thermodynam- 
ics, chemotaxis, gravity, long-range interactions 



1 Introduction 

In biology, many microscopic organisms (bacteria, amoebae, endothelial cells,...) or even so- 
cial insects (like ants) interact through the phenomenon of chemotaxis pQ. These organisms 
deposit a chemical (pheromone, smell, food,...) that has an attractive effect on the organisms 
themselves. Therefore, in addition to their diffusive motion, they move preferentially along the 

1 The case of repulsive chemotaxis due to a "poison" can also be of interest and will be considered in a future 
contribution. 
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gradient of concentration of the chemical they secrete (chemotactic flux). When chemotactic 
attraction prevails over diffusion, this process can lead to a "chemotactic collapse" (see [2] for 
a review) resulting in the aggregation of the organisms. In this way, some structures can form 
like clusters (clumps) or even network patterns (filaments). Therefore, the chemotactic interac- 
tion can explain several features of the morphogenesis of biological colonies. The chemotactic 
aggregation of biological populations is usually described in terms of the Keller-Segel model 
[3j. This is a parabolic model consisting in two coupled differential equations. The first equa- 
tion is a drift-diffusion equation describing the evolution of the concentration of cells and the 
second equation is a reaction-diffusion equation with terms of source and degradation describ- 
ing the evolution of the concentration of the secreted chemical. This model ignores inertial 
effects and assumes that the drift velocity of the organisms is directly induced by a chemotactic 
"force" proportional to the concentration gradient of the chemical. The Keller-Segel model 
can reproduce the formation of clusters (clumps) by chemotactic collapse [4-19]. This reflects 
experiments on bacteria like Escherichia coli or amoebae like Dictyostelium discoi'deum exhibit- 
ing pointwise concentration [3J. However, parabolic models fail at describing the formation of 
network patterns (filaments). These filaments are observed in experiments of capillary blood 
vessels formation [2U|. They correspond to the spontaneous self-organization of endothelial cells 
during vasculogenesis, a process occuring during embryologic development. In order to account 
for these structures, more general models of chemotaxis have been introduced [21], [221 [22] . They 
have the form of hydrodynamic (hyperbolic) models taking into account inertial effects. These 
models can reproduce the formation of filaments that are interpreted as the beginning of a 
vasculature. This phenomenon is responsible of angiogenesis, a major factor for the growth of 
tumors [24] . Interestingly, these filaments share some analogies with the large-scale structures 
in the universe that are described by similar hydrodynamic equations [25l [26] . 

Recently, we have introduced a general kinetic model of chemotactic aggregation based on 
generalized stochastic processes, non linear mean field Fokker-Planck equations and generalized 
thermodynamics [23] . From these kinetic equations, we have derived a hydrodynamic model of 
the form 

(1) ! + v-(H = o, 

(2) ^ + ( u . V )u=--Vp + Vc-£u, 

at p 

3c 

(3) — = D c Ac -kc + hp. 

It involves a pressure force — Vp, where p = p(p) is a barotropic equation of state that can 
take into account several effects like anomalous diffusion or the fact that the particles do 
not interpenetrate. It also involves a friction force — £u which measures the importance of 
inertial effects. For £ = 0, we recover the hyperbolic model introduced by Gamba et al. [21] . 
For £ — > +oo, we can neglect the inertial term in the momentum equation (T5]) leading to 
£u ~ — ^Vp + Vc (overdamped limit). Substituting this relation in the equation of continuity 
(JTJ) , we obtain the generalized Keller-Segel model 

(4) ^ = V ■ [ X (Vp - pVc)] , 

3c 

(5) — = D c Ac -kc + hp, 
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where x — V£- Interestingly, this model of chemotaxis is similar to a model of self-gravitating 
Langevin particles [27] described by the damped Euler-Poisson system 

(6) ^ + V-(pu) = 0, 



(7) 



du 



+ (u- V)u 



1 

P 



-Vp - V$ - £u, 



(8) 



A$ = S d Gp. 



For £ = 0, it reduces to the barotropic Euler-Poisson system [28] and for £ 
the generalized Smoluchowski-Poisson system 



-oo, we obtain 



(9) 



at 



- (Vp + pV$) 



(10) A$ = S d Gp. 

In this analogy, we see that the concentration — c(r, t) of the chemical plays the same role 
as the gravitational potential $(r, t). In biology, the interaction is mediated by a material 
substance (the secreted chemical) while the physical interpretation of the gravitational potential 
in astrophysics is more abstract q The hydrodynamic equations ([I])-© or fl6^- (fT0l) involving a 
barotropic equation of state, a long-range potential of interaction and a friction force, have been 
introduced by Chavanis j29j EQ] at a general level. It was indicated that they could provide 
generalized models of chemotaxis and self-gravitating Brownian particles. 

The main difference between the chemotactic model (DO)-© and the gravitational model 
(JS])- (fIU]) concerns the field equations <^ and <$E§. In astrophysics, the gravitational potential 
is determined instantaneously from the density of particles through the Poisson equation (jSJ). 
In biology, the equation ([3]) determining the evolution of the chemical is more complex and 
involves memory terms. The chemical diffuses with a diffusion coefficient D c , is produced by the 
organisms at a rate h and is degraded at a rate k. Because of the term dc/dt, the concentration 
of the chemical at time t depends on the concentration of the organisms at earlier times. In 
this paper, we shall consider simplified models where the term dc/dt can be neglected. This 
is valid in a limit of large diffusivity of the chemical D c — > +oo [6]. We first consider the case 
where there is no degradation of the chemical (k = 0). Then, assuming h = XD C and taking 
the limit D c — > +oo with A = 0(1), one gets (see Appendix C of [23] ) 

(11) Ac=-\(p-p), 

where p = {1/V) J pdr = M/V is the average value of the density which is a conserved quantity. 
In that case, the concentration of the chemical is given by a Poisson equation which incorporates 
a sort of "neutralizing background" (played by p) like in the Jellium model of plasma physics 
[31] . Note that a similar term also arises in cosmology when we take into account the expansion 
of the universe and work in a comoving system of coordinates [25]. We shall thus refer to this 
model as the "Newtonian model". Then, we consider the case of a finite degradation rate. 

2 The notion of "force at distance" in the Newtonian theory has been criticized at several occasions in the 
history of physics and replaced by the notion of curved space-time in the Einsteinian theory. 
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Assuming h = XD C , k = k\D c and taking the limit D c — > +00 with A = 0(1) and k = 0(1), 
one gets (see Appendix C of [23] ) 




If we take formally k$ = 0, we obtain a Poisson equation similar to Eq. ([8]) where — c(r, t) 
plays the role of <3>(r, t) and A plays the role of the gravitational constant SdG (we recall that 
the geometrical factor Sd is the surface of a unit sphere in d dimensions). However, Eq. (IT2T) 
has been derived for fc 7^ (for /c — we get Eq. (jlip ). This implies that the interaction is 
shielded on a typical distance A;^" 1 . This is similar to the Debye shielding in plasma physics, 
to the Rossby shielding in geophysical flows or to the Yukawa shielding in nuclear physics. We 
shall refer to this model as the "Yukawa model" . 

In this paper, we perform a detailed linear dynamical stability analysis of the chemotactic 
model ([I])-©. Specifically, we study the stability of an infinite and homogeneous distribution 
of cells against chemotactic collapse. This is similar to the classical Jeans stability analysis 
for the barotropic Euler- Poisson system [28]. Indeed, the "chemotactic collapse" of biological 
populations is similar to the "gravitational collapse" in astrophysics (Jeans instability). There 
are, however, two main differences with the classical Jeans analysis. The first difference is the 
presence of a friction force — £u in the Euler equation. As we shall see, this does not change the 
onset of the instability but this affects the evolution of the perturbation. The second difference 
arises from the different nature of the field equations ([3]) and (181) . We recall that, in gravitational 
dynamics, an infinite and homogeneous distribution of matter with p = est and u = is not 
a stationary solution of the barotropic Euler- Poisson system flS])-© because we cannot satisfy 
simultaneously the condition of hydrostatic equilibrium Vp(p) + pV$ = reducing to = 
and the Poisson equation A$ = S^Gp 7^ 0. This leads to an inconsistency in the mathematical 
analysis when studying the linear dynamical stability of such a distribution: this is called the 
"Jeans swindle" [28] Q By contrast, there is no "Jeans swindle" in the chemotactic problem! 
Indeed, an infinite and homogeneous distribution of cells is a steady state of the equations of 
motion ([I])-© corresponding to the condition kc = hp. For the "Newtonian model" (ITT]) , this 
condition becomes p = ~p and for the "Yukawa model" (IT2l . it becomes k\c = Xp. 

In this paper, we study in detail the onset of the "chemotactic instability" and its develop- 
ment in the linear regime. This study was initiated in [31] at a general level, i.e. taking into 
account the term dc/dt in Eq. ([3]) and allowing the coefficients in Eqs. ([I])-® to depend on 
the concentration. However, this study focused on the unstable modes and did not analyze in 
detail the evolution of the stable modes. In the present paper, we make a complete study of 
both stable and unstable modes but we restrict ourselves to the simplified models (II ip and (j!2p . 
In the "Newtonian model" (ITT]) , the only difference with the Jeans analysis is the presence of 
the friction force £. In the "Yukawa model" ()12p . the differences with the Jeans analysis are 
due to the effects of the friction £ and of the shielding length k^ 1 generated by the degradation 

3 One possibility to avoid the Jeans swindle is to study the linear dynamical stability of an inhomogeneous 
distribution of matter in a finite domain (box) |32j . Alternatively, in cosmology, the "Jeans swindle" is cured 
by the expansion of the universe [25j . Indeed, if we work in a comoving system of coordinates, the usual Poisson 
equation A$ = AnGp is replaced by an equation of the form A(f> = 4irGa(t) 2 [p(x,t) — pb{t)\ where the density 
p(x, t) is replaced by the deviation p(x, t) — pb(t) to the mean density |25| . Then, an infinite and homogeneous 
distribution of matter with p = p\j and <p = is a steady state of the equations of motion from which we can 
develop a rigorous stability analysis. The expansion of the universe introduces a sort of neutralizing background 
in the Poisson equation. Interestingly, the same effect arises in the chemotactic model (|TTj) for a completely 
different reason. Note finally that, in early models of cosmology, some authors including Einstein himself have 
modified the gravitational Poisson equation to the form A<E> — A<& = AirGp by including a shielding term |33j . 
This transformation was done in order to obtain a static homogenous and isotropic universe. As we have seen, 
a similar shielding effect arises naturally in the chemotactic model ([T2|) due to the degradation of the chemical. 




(12) 



Ac — k?.c 



■ c - 



Xp. 
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of the chemical. We show that the system is always stable for 

/ Xp\ 

(13) c s > (c s ) cr it = y~^2 J ' 

where c s = {dp/ dp) 1 / 2 is the "velocity of sound" in the medium (for specific equations of state, 
discussed in Sec. HI we can express the stability criterion ffl3|) in terms of the density of cells). 
Therefore, the system is stable if the velocity of sound is above a certain threshold fixed by the 
shielding length k . By contrast, for c s < (c s ) cr i t , the system is unstable for wavenumbers 

(14) k < k m EE \J k 2 j - kl 

where kj = (Xp/c 2 ) 1 / 2 is the Jeans wavenumber. In the Newtonian model, the condition k = 
implies (c s ) cr i t = +oo, so that the system is always unstable to perturbations with sufficiently 
large wavelengths k < kj. These results are independent on £. The friction term only affects the 
evolution of the perturbation. For k < k m (ko), the perturbation grows exponentially rapidly, 
for k m (ko) < k < k c (£,,ko) (where k c is a friction-dependent wavenumber defined in the text) 
it is damped exponentially rapidly without oscillating and for k > fc c (£, fco) it presents damped 
oscillations. More precisely, we determine the growth rate of the unstable modes and the 
damping rate, and oscillation frequency, of the stable modes as a function of £ and k^ 1 . Owing 
to the above mentioned analogy between chemotaxis and gravity, our stability analysis also 
applies to self-gravitating Langevin particles [27] provided that we make the "Jeans swindle" . 



2 Jeans-type instability for a Newtonian potential 

In this section, we study the linear dynamical stability of an infinite and homogeneous stationary 
solution of the fluid equations 

(15) ^ + V • (pu) = 0, 

(16) -^ + (u-V)u = — Vp + Vc-£u, 

at p 

(17) Ac=-X(p-p). 

We consider an infinite and homogeneous distribution of cells p = ~p, with no velocity u = 
and no chemical c = 0. This is an exact stationary solution of the fluid equations ( fTBl -lfTTj). 
Linearizing Eqs. ( fT5i) -( fiTj) around this steady state and writing the perturbation in the form 
Sf(r,t) ~ e *(k-r-w^ we reac iiiy obtain the dispersion relation [27] : 

(18) u(u + z£) = c 2 k 2 - Xp, 

where we have introduced the equivalent of the velocity of sound c 2 = p'(~p). Setting a = —iu, 
so that 5f oc e at , the dispersion relation can be rewritten 

(19) a 2 + £a + c 2 s k 2 - Xp = 0. 

The solutions are o± = (-£ ± y/~K)/2 with A(fc) = £ 2 - A(c 2 k 2 - Xp). If c 2 k 2 - Xp < 0, then 
A > £ 2 > and the system is unstable since a + = (— £ + vA)/2 > 0. If c 2 k 2 — Xp > 0, either 
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(i) A < implying R e (a) = -£/2 or 
Therefore, the system is unstable if 

(20) k < 



< A < £ 2 , implying a± < 0, so the system is stable. 



Xp 



1/2 



= k A 



and stable otherwise. The critical value kj is similar to the Jeans wavenumber in astrophysics. 
We note that the threshold of instability does not depend on the friction parameter £. Note 
also that for negative chemotaxis (chemorepulsion) obtained by replacing +Vc by — Vc in Eq. 
(I16p . an infinite and homogeneous distribution of particles is always stable. 

If we consider the case £ = (Euler), the fluid equations ([T5"|) - (fT7")) are similar to the Euler- 
Poisson system and the dispersion relation becomes 



(21) 



For k > kj, the perturbation undergoes undamped oscillations with pulsation to = c s (k 2 —k 2 ) 1 ^ 2 . 
For k < kj, the perturbation increases exponentially rapidly with a growth rate 7 = c s (k 2 — 
k 2 ) 1 / 2 . For c s = 0, kj — > +00 and the system is unstable for all wavenumbers. The growth rate 
of the perturbation is 7 = \fXp independent on k. For c s — > +00, kj — >• and the system is 
stable for all wavenumbers. The pulsation is u = c s k. If we now consider the case £ — > +00 
(Smoluchowski), the fluid equations (fl5j) - (fTTj) reduce to the generalized Smoluchowski-Poisson 
system 



(22) 



dt 



^(Vp - pVc) 



(23) Ac=-X(p-p), 
and the dispersion relation becomes 

(24) i^uj = c 2 s k 2 - Xp. 

For k > kj, the perturbation decays exponentially rapidly with a damping rate 7 = —c 2 {k 2 — 
k 2 )/^. For k < kj, the perturbation increases exponentially rapidly with a growth rate 7 = 
c 2 (k 2 — k 2 )/^. For c s = 0, kj — > +00 and the system is unstable for all wavenumbers. The 
growth rate of the perturbation is 7 = Ap/£ independent on k. For c s — > +00, kj — > and the 
system is stable for all wavenumbers. The damping rate is 7 = —c 2 k 2 /^. 

Let us now consider the general case of an arbitrary friction. There are two relevant 
wavenumbers in the problem: the Jeans wavenumber (12"U1) and the wavenumber 

(25) k c =(k 2 +k 2 ) 1/2 , 
where 

(26) h . JL, 

is a wavenumber constructed with the friction coefficient and the velocity of sound. We can 
define a dimensionless number 

k d \ 2 e 



(27) F=[f j 



AXp 



6 



2 




12 3 4 

k 2 

Figure 1: In summary, a homogeneous distribution is unstable for k < kj and stable for k > kj. 
For k < kj, the perturbation grows exponentially rapidly. For kj < k < k c , the perturbation 
is damped exponentially rapidly without oscillating. For k > k c , the perturbation undergoes 
damped oscillations. We have taken F — 1, kj — 1 and £ = 2. 

which measures the strength of the friction force (a similar parameter was introduced in [27] 
for inhomogeneous distributions). It is independent on the equation of state p(p) and it can 
be written F ~ (t,t D ) 2 where ~ l/^/^pX is a typical dynamical time. Thus, y/F is the 
ratio of the dynamical time on the friction time r ~ l/£. In terms of this parameter, the 
wavenumber (1251) can be written k c = kj\J\ + F. The behaviour of the perturbation can be 
analyzed in terms of these wavenumbers: (i) If A < 0, the perturbation undergoes damped 
oscillations with pulsation to = c s (k 2 — k 2 ) 1 ^ 2 and decay rate 7 = — £/2. This stable regime 
corresponds to wavenumbers k > k c . (ii) If < A < £ 2 , the perturbation decays exponentially 
rapidly with a damping rate 7 = — £/2 + c s (k 2 — k 2 ) 1 ^ 2 without oscillating. This stable regime 
corresponds to wavenumbers kj < k < k c . For k = kj, we have 7 = and for k = k c , we 
have 7 = — £/2. (iii) If A > £ 2 , the perturbation increases exponentially rapidly with a growth 
rate 7 = -£/2 + c s (k 2 c - k 2 ) 1 ' 2 . This unstable regime corresponds to wavenumbers k < kj. 
The growth rate is maximum for k* = and its value is 7* = — £/2 + c s k c . These results are 
summarized in Fig. [TJ 

It is interesting to determine how the results depend on the friction parameter and on the 
velocity of sound. To simplify the notations, we define T = c 2 . Then, we obtain 

(28) kj(T)=(^y\ 

(29) k c (T,0 = kj(T)VT+F, 



(30) fc, = 0, | 7 ,(e) = -l + yi + i 

1/2 

(k < k c ). 



(31) 



1 

F 



7 



2 




(32) 



>4 



1/2 



7 



-e/2, (A; > k c ). 



For T = 0, we find that kj — > +oo so that the system is unstable for all wavenumbers. The 
growth rate of the perturbation is 



(33) 



27 



-l 



independent on fc. For T — > +oo, if £ is finite, we get kj = k c = so that the system is stable 
for all wavenumbers. The pulsation is u = VTk and the damping rate 7 = — £/2. 
For £ = 0, we find that k c = kj and 



(34) 



7(fc,T) = 



1- 



fc 



1/2 



(A; < fcj) 



(35) 



cu(fc, T) = a/Xp 



fc 

^7 



1/2 



7 = 0, 



{k > kj). 



These results are summarized in Fig. [21 For T = 0, we find that kj — > +00 so that the system is 
unstable for all wavenumbers. The growth rate of the perturbation is 7 = \/Xp. For T — > +00, 
we get /cj = so that the system is stable for all wavenumbers. The pulsation is u — y/Tk. 
For £ — > +00, we find that A; c — > +00 and 



(36) 



T (fc,T) 



1 - 



k 2 j 



These results are summarized in Fig. [3j For T = 0, we find that kj — > +00 so that the 
system is unstable for all wavenumbers. The growth rate of the perturbation is 7 = Ap/£. For 
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Figure 3: The limit £ — > +00. We have taken kj — 1 and Xp = 1 

T — > +00, we get fcj = so that the system is stable for all wavenumbers. The damping rate 
is 7 = -Tk 2 /i. 

If we now consider the stability problem of Eqs. (fT5l) - (|T7|) in a two-dimensional periodic 
domain of size L, the wavenumbers can be written k = x( m ' n ) wriere m, n are positive integers 
with (m,n) 7^ (0,0). In that case, the condition of instability (f2"0"|) becomes 

(37) m 2 + n 2 < ApL 2 



4:7T 2 C 2 S ' 

A necessary condition of instability is therefore XpL 2 /(4ti 2 c 2 ) > 1. For an equation of state of 
the form p(p) = pT, where T plays the role of a temperature, the velocity of sound is c 2 = T 
and the necessary condition of instability can be written 

(38) T < T c = 

where T c is a critical temperature. For T > T c , there is no chemotactic collapse: the "gas" of 
cells remains spatially uniform and diffuse. For T < T c , the distribution of cells is unstable and 
the number of unstable modes increases as T decreases yielding more and more clusters. This 
instability has been illustrated numerically in [23] by solving the iV-body equations of motion 
in a two-dimensional periodic domain. 

3 Jeans-type instability criterion for a Yukawa potential 

We now consider the linear dynamical stability of an infinite and homogeneous stationary 
solution of the fluid equations 

(39) ^ + V • (pu) = 0, 

(40) -7^- + (u ■ V)u = — Vp + Vc - f u, 

ut p 
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(41) Ac - k 2 c= -Xp. 

Comparing with Eq. (ITT]) , we see that the Laplacian is replaced by the operator A — k^. This 
implies that the interaction, mediated by the concentration c of the chemical, is screened on a 
distance k^ 1 where k = (k/D^ 1 / 2 [k should not be confused here with the wavenumber]. A 
uniform distribution of cells and secreted chemicals whose concentrations satisfy the relation 
k\c = Xp is an exact stationary solution of Eqs. (l39l) - fj4~Tl) . Considering a small perturbation 
around this steady state, we find that the dispersion relation replacing Eq. ffT9~]) is 



The solutions are a± = (— £ ± v A)/2 with 
(43) A(k) =e- 4^ 



Xp 



k 2 + kl 

Repeating the arguments following Eq. (fl9|) . we find that the system is unstable if 

(44) c 2 < X ? 

{ ' s k 2 + k 2 ' 

and stable otherwise. A necessary condition for instability is that 

(45) c 2 s < (c 2 ) crU ee 

When this condition is fulfilled, the range of unstable wavelengths is 



(46) k < k m = y/k* - k 2 . 

We see that for ko ^ 0, the instability is shifted to larger wavelengths than the Jeans length. 
If we consider the case £ = (Euler), the dispersion relation becomes 



For k > k m , the perturbation undergoes undamped oscillations with pulsation uj = \/ — a 2 . For 
k < k m , the perturbation increases exponentially rapidly with a growth rate 7 = y/a 2 . If we 
consider the case £ — > +00 (Smoluchowski), the dispersion relation becomes 

k 2 - P 

(48) {„ = <ik' 2 ^J. 

For > fc m , the perturbation decays exponentially rapidly with a damping rate 7 = o < 0. 
For < fc m , the perturbation increases exponentially rapidly with a growth rate 7 = cr > 0. 

Let us now consider the general case of an arbitrary friction. We introduce a critical 
wavenumber k c defined by 



(49) 2k 2 c = k 2 m + k 2 + ^(k 2 m + k 2 ) 2 + Ak 2 k% 

This expression generalizes Eq. (J2"5]) to the case k ^ 0. We also introduce the wavenumber k s 
defined by 



(50) 2k 2 s = -{kl + k 2 d ) + sj{k 2 m + k 2 ) 2 + 4k 2 k 2 d . 
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0.5 



1.5 



Figure 4: Growth rate and pulsation as a function of the wavenumber. We have taken T/T c = 
1/2, F = 0.05, k = 1 and f = 2 

The behaviour of the perturbation can be analyzed in terms of these wavenumbers: (i) If A < 0, 
the perturbation undergoes damped oscillations with pulsation 



(51) 

or equivalently 
(52) 



£ jk 2 (k 2 -k 2 n 
2 



to 



(fc 2 -fc e 2 )(fc 2 + fc 2 ) 
k 2 + k 2 



1/2 



and decay rate 7 = — £/2. This stable regime corresponds to wavenumbers k > k c . (ii) If 
< A < £ 2 , the perturbation decays exponentially rapidly with a damping rate 



(53) 

or equivalently 
(54) 



k 2 (k 2 -k 2 J 
k 2 (k 2 + k 2 ) 



7 



i +Cs 



(fc e 2 -fc 2 )(fc 2 + fc 2 ) 

k 2 + k 2 



1/2 



without oscillating. This stable regime corresponds to wavenumbers k m < k < k c . For k = k c , 
we have 7 = — £/2 and for k = k m we have 7 = 0. (iii) If A > £ 2 , the perturbation increases 
exponentially rapidly with a growth rate 7 given by Eq. fl53|) . This unstable regime corresponds 
to wavenumbers k < k m . The growth rate is maximum for 



(55) 

and its value is 
(56) 



1* 



K = \Jhikj - koj, 



c s Jk 2 + (kj-k y 
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In summary, a homogeneous distribution is unstable for k < k m and stable for k > k m . For 
k < k m , the perturbation grows exponentially rapidly For k m < k < k c , the perturbation 
is damped exponentially rapidly without oscillating. For k > k c , the perturbation undergoes 
damped oscillations. These results are summarized in Fig. HI 

It is interesting to determine how the results depend on the friction parameter and on the 
velocity of sound (see Figs. and [5]) • To simplify the notations, we define T = c 2 . Then, we 
obtain 

/ T \l/2 /rp \i/2 /T\ 1/2 

(57) kj(T)/k = , k m (T)/k = ^ - lj , k d (£, T)/k — VF \ t^J , 



(58) 



2fc c 2 (e,T)A 2 = (l + F)^-l 



n 2 



T 
T 



(59) 



2k 2 (^T)/k 2 = 1-(1 + F)^ + 



T 
T 



(60) 



K(T)/k 



1/2 



1/2 



(61) 



- 7 *(e,T) = -i + 



1 

1+ F 



1/2 



(62) 



2 „ . / k 2 [k 2 - k 2 (T c /T - 1)1 



fc 2 F(T c /T)(^ + fc 2 



(A; > A; c 



(63) 



e 7(,e ' } " + V kiF(r c /Tw + ki)' {k<kc 



For T = T r , we have 



(64) 



fcj = k , k m = 0, fc d = k \/F, 



(65) 



2fc 2 /fc 2 = F + VF 2 + 4F, 2A; 2 /A; 2 = — F + VF 2 + 4F, 



(66) 



2 

f V W^ 2 + kl) 



1, ^ = _ 2' ( k>kc ">> 



(67) 



:7 = -1 + a/1 



fc 4 



fc 2 F(£; 2 + fc 2 )' 



(A; < fc c ). 
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Figure 5: Dependence of the characteristic scales with the temperature. The growth rate 
is maximum for k = k*(T) and the system is stable for k > k m (T). Oscillations appear 
for k > k c (T,F). For T — > 0, we have k 2 Jk\ ~ (1 + F)T C /T and for T — > +oo, we have 
fc^/fcg ~ FT C /T. Note that for F = 0, k c = k m for T < T c and fc c = for T > T c . 
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Figure 6: Maximum growth rate as a function of the temperature for different values of the 
friction parameter F. 
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The system is stable for all wavelengths. The wavenumber marking the appearance of oscilla- 
tions is k c (F). It behaves like k c /k ~ F 1 / 4 for F -> and k c /k ~ F 1 ' 2 for F -> +00. 

For T = 0, we find that fc m — > +00 so that the system is always unstable. The growth rate 
of the perturbation is given by 



(68) 



2 T- 



1 + 



k 2 



F(k 2 + k 2 ) 

It increases and tends asymptotically to its maximum value 



(69) 



For T > T c , the system is always stable (k 2 2 < 0). For k > k c , the perturbation undergoes 
damped oscillations with pulsation floTl) and decay rate 7 = — £/2. For fc < fc c , the perturbation 
decays exponentially with rate (153]) tending to 7 = for k = 0. For T — > +00 and F finite, 



(70) 



0. 



0. 



0. 



The perturbation oscillates with a pulsation u = VTk, and is damped with a rate 7 
The case F — >• +00 is treated in Sec. 



3.1 The case ^ = 

For £ = (Euler) and T < T c , we have 



(71) 



- T \ 1/2 
MT) = fcc(r) = fco ( - 1 



= fc. = 0. 



For k < k m , the system is unstable and the growth rate is 



(72) 

It is maximum for 
(73) 

with value 
(74) 



7 = VT 



k\k 2 m - g 
fc 2 + k 2 



1/2 



k*(T)/k 



1/2 



1/2 



7*(T) = /cqv 7 ^ 



1/2" 



For k > k m , the system is stable and the perturbation undergoes oscillations with pulsation 



v k\k 2 -k 2 m ) 



(75) 



1/2 



k 2 + k 2 

These results are summarized in Fig. [71 For T = T c , k m = 0. The system is always stable and 
the pulsation is 



(76) 
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Figure 7: Growth rate and pulsation as a function of the wavenumber for £ = 0. We have taken 
T/T c = 1/2, Ap = l, fc = l. 

For T = 0, fc m ~ ko(T c /Ty/ 2 — > +oo. The system is always unstable and the growth rate is 

kok 



(77) 



1=VT< 



Vk^Tkl' 



It increases and tends, for k — > +oo, to its maximum value 7* = vTc&o = y\p. For £ = and 
T > T r , we have 



(78) 



7; 

T 



k d = k c , = 0. 



The system is always stable and the pulsation is 

'k 2 (k 2 + k 2 s y 1/2 



(79) uj = VT 

For T — > +00, k s — > fco and ~ vTA;. 



fc 2 + fcg 



3.2 The case £ -> +00 

For £ — »• +00 (Smoluchowski) and T < T c , we have 



(80) 



k m {T) = fc ( ^ - 1 



1/2 



^ c ~ kd — »■ +00, fc s = 0. 



The rate of the exponential evolution is 



(81) 
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r fc 2 (*&-fc s 

£ /c 2 + fcg 
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Figure 8: Growth rate as a function of the wavenumber for £ — > +oo. We have taken T/T c = 
1/2, \p=l, k = l. 

For k > k m , the perturbation is damped (stable) and the decay rate behaves like 7 ~ — (T/^)k 2 
for k — > +00. For k < k m the perturbation increases (unstable). The growth rate is maximum 
for 



(82) 

with value 

(83) 



K(T)/k 



1/2 



1/2 



7*(T) 



T 



1/2 



These results are summarized in Fig. [HI 



For T = T c , k m = 0. The system is always stable and the perturbation decreases with a 



decay rate 
(84) 



T r k 4 



7 



e ^ + ki 



For T = 0, k m —> +00. The system is always unstable and the growth rate is 



(85) 



7 



k(yT c k 2 
£ k 2 + k 2 



It increases and tends, for k —> +00, to its maximum value 7* = k 2 T c /^ = Ap/£. For £ — > +00 
and T > T r , we have 



(86) 



k 2 m = k 2 ( 5 - 1 



T 



k c ~ k d ^ +00, 



k s = 0. 



The system is always stable and the decay rate is 
(87) ^_Tk 2 (k 2 m -k 2 ) 

For T - 



i k 2 + k 2 



< 0. 



-00, k 2 m 



-&q and 7 = — jk 2 
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4 Particular equations of state 



The hydro dynamic equations ([I])-© or the generalized Keller-Segel model (jl])-© incorporate 
a pressure force — Vp(p) associated with a barotropic equation of state p(p). At equilibrium, 
the system satisfies a relation of the form 

(88) Vp = pVc. 

This corresponds to a condition of hydrostatic balance between the pressure force and the 
chemotactic attraction. For inhomogeneous systems at equilibrium, the density is a function 
p = p(c) of the concentration of the chemical obtained by integrating Eq. (JSSl) . We have the 
identities 

/ x /">'(» , P'(p) 1 

89 / ^dx = c, ZW_ = —, p'(c) = p. 

J x p p'(c) 

Since p'(p) > in ordinary circumstances, this implies that the density is an increasing function 
of the concentration of the chemical, i.e. p'(c) > 0. On the other hand, for homogeneous 
systems, the condition of stability (fT3l) can be written 

(«» <-f > A 

P K o 

This relation determines the range of densities p for which the system is stable, depending on 
the form of the equation of state p(p). In particular, the system is stable whatever the value of 
the density if mm p [p'(p)/p\ > X/k 2 , and unstable (for sufficiently large wavelengths) whatever 
the value of the density if max p [p'(p)/p] < A/A;q. The pressure force in Eq. can take into 
account different effects such as anomalous diffusion or close packing effects. Typically, three 
kinds of pressure law p(p) have been considered in the chemotactic literature: 
(i) In the standard case, the pressure is a linear function of the density 

(91) p = pT. 

This is similar to an isothermal equation of state where T is an effective temperature [H] . When 
this law is substituted in the drift-diffusion equation (jlj) we recover the standard Keller-Segel 
model [3]: 

(92) ^ = V .(DVp-xpVc), 

where we have introduced the diffusion coefficient D through the Einstein relation D = yT . 
The equilibrium state is the Boltzmann distribution p = Ae^ c . For the pressure law (|9T|) . the 
velocity of sound has the constant value c 2 s = T. Therefore, the stability criterion f|T3|) can be 
rewritten in the form 

(93) t>T c =^§. 

For a fixed density p, it defines a critical temperature T c below which the system is unstable. 
Alternatively, for a fixed temperature T, the stability criterion fTTBl can be rewritten in the 
form 

k 2 T 

(94) p < p m t = 
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The system becomes unstable above a certain critical density p cr u- 

(ii) In [15], [35], we have proposed to take into account anomalous diffusion by using an 
equation of state of the form 

(95) p = Kp 1 . 

This is similar to a polytropic equation of state where K plays the role of a polytropic tempera- 
ture. When this law is substituted in the drift-diffusion equation (jlj) we obtain the generalized 
Keller-Segel model studied in [TBI 135] : 

(96) ^ = V ■ [ X (KVp^ - pVc)] . 

The equilibrium state is the Tsallis distribution p = [A + (7 — \)/{K , y)(^'^ f ~ x ' [36]. For the 
pressure law (1931) . the velocity of sound has the value c 2 s = K^p 1 ' 1 . Therefore, the stability 
criterion ([TBI can be rewritten in the form 

(97) K > K c = 



For a fixed density p, it defines a critical polytropic temperature K c below which the system 
is unstable. Alternatively, for a fixed polytropic temperature K, we can express the stability 
criterion (TIBl) as a function of the density. We need to distinguish three cases (see Figs. [9] and 
[TO]) : (a) For 7 < 2, the stability criterion can be written in the form 



(98) p<p 



crit 



2-7 



The system becomes unstable above a critical density, (b) For 7 > 2, the stability criterion can 
be written in the form 



(99) p>p 



crit 



A A -<- 2 



The system becomes unstable below a critical density, (c) For 7 = 2, the stability criterion can 
be written in the form 

(100) K > K c = A 

The instability threshold is independent on the density. 

(iii) As a result of chemotactic collapse, the standard Keller-Segel model can lead to finite 
time singularities and Dirac peaks [21 HH EE]. In reality, these singularities are unphysical 
because the cells have a finite size and cannot be compressed indefinitely. Therefore, in more 
realistic models, we expect that the pressure p(p) tends to zero for low densities p — > and 
rapidly increases for large densities. This takes into account the fact that cells do not interpen- 
etrate due to their finite size and this prevents overcrowding. In [37] , one of us has proposed to 
take into account volume filling and finite size effects by using an equation of state of the form 

(101) p=-Ta \n(l-p/a ). 

For low densities p — > 0, we recover the linear equation of state p = pT and for high densities, 
close to the maximum allowable density ctq, the pressure rapidly increases and diverges when 
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p 



Figure 9: Graphical construction determining the instability threshold (expressed in terms of 
the density) for a polytropic equation of state with index 7. The system is stable if c 2 s = 
K'jp 1 ' 1 > \pjk\ and unstable to sufficiently large wavelengths otherwise. We have taken 
K — 1 and A/fcg = 1. 




Figure 10: Stability diagram of a spatially uniform polytropic gas with an attractive Yukawa 
potential of interaction. The curves represent K c (p) or p crit (K) and separate the stable region 
(upper region) from the unstable one (lower region). For 7 = 1, we recover the isothermal case 
with K = T. We have taken A/ifeg = 1. 



19 




Figure 11: Graphical construction determining the instability threshold (expressed in terms of 
the density) for an equation of state of the form ( 110 II) with temperature T. The system is 



stable if c 2 s = T/(l 
We have taken ctq = 



- p/<7q) > Ap/A;Q and unstable to sufficiently large wavelengths otherwise. 
1 and = 1 leading to a critical temperature T* = 1/4. 



0.4 



0.3 



H0.2 



Stable 




Figure 12: Stability diagram of a spatially uniform gas described by the pressure law (11011) 
with an attractive Yukawa potential of interaction. The curve represents T c (p) or p±(T) and 
separates the stable region from the unstable one. There exists a critical point T* above which 
the homogeneous phase is always stable whatever the density. Below T*, the homogenous phase 
is stable for p < p_, unstable for p_ < p < p + and stable again for p > p + . This corresponds 
to a reentrant phase. We have taken er = 1 and A/fcg = 1 leading to a critical temperature 
T* = 1/4. 
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p — > a . If a represents the typical size of the cells, we have a ~ l/a d where d is the dimension of 
space. When this law is substituted in the drift-diffusion equation (j4j) we obtain the generalized 
Keller-Segel model studied in 



The steady state of this equation is a Fermi-Dirac distribution in physical space p = cr /(l + 
Ae _/3c ) putting an upper bound on the density: p < a 0. For the pressure law (jlOip . the 
velocity of sound is c 2 s = T/(l — p/cr ). Therefore, the stability criterion (fTBl can be rewritten 
in the form 

(103) T >T C = 

For a fixed density p, it defines a critical temperature T c below which the system is unstable. 
Alternatively, for a fixed temperature T, we can express the stability criterion ( TTBl) as a function 
of the density. The system is stable if 

(104) p 2 -a p+^>0, 

and unstable to large wavelengths otherwise. The discriminant of this equation is A = o-q — 
4T<7 ^o/^- If A < 0, corresponding to 



(105) T > r„ 



Aero 




the system is stable whatever the value of the density (see Figs. [TT1 and fT2l) . Alternatively, if 
T < T*, the system is stable for p < p_ and p + < p < o" and unstable for p_ < p < p + where 

(106) p ± = 

For low densities (p o" ), the homogenous phase is stable because the chemotactic attraction 
is not strong enough to overcome diffusive effects (like for a low density isothermal gas (i)) and 
for high densities (p — > er )> the homogeneous phase is stabilized by pressure effects due to close 
packing. Other generalized Keller-Segel models of chemotaxis are discussed in |39| in relation 
with nonlinear mean field Fokker-Planck equations. 



5 Conclusion 

In this paper, we have studied the chemotactic instability of an infinite and homogeneous 
distribution of cells whose dynamics is described by the hydrodynamic equations ([I])-©. We 
have shown the analogy with the classical Jeans instability in astrophysics. This close analogy 
between two systems of a very different nature (stars and bacteria) is very intriguing and 

4 Equation (|102[) can be viewed as a generalized mean field Fokker-Planck equation [52) with a constant 
mobility and a nonlinear diffusion. A related model, corresponding to a constant diffusion and a variable 
mobility \(p) = xO-~P/ cr o) vanishing above the close packing value o"o, has been considered in [38l [29l [37] . The 
two models have the same steady states and are associated with the same free energy. They present therefore 
the same general properties. The details of the evolution may, however, be different in the two models. 
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deserves to be developed and emphasized. As is well-known, an infinite and homogeneous 
distribution of stars is not a stationary state of the gravitational Euler-Poisson system [28J. 
However, if we make the Jeans swindle (which is made in any textbook of astrophysics), the 
equations for the linear perturbations are the same as in the biological problem ([I])-© when 
£ = ko = 0. Therefore, the two systems are really analogous. The main differences between the 
chemotactic problem and the Jeans problem are due to the presence of (i) a friction force — £u 
in the Euler equation (J2J) and (ii) a shielding length k$ l in the equation (j3J) determining the 
potential of interaction (played here by the concentration of the secreted chemical). We have 
studied the effect of these terms in detail. This leads to a generalization of the Jeans instability 
analysis. The shielding length determines a critical velocity of sound (c s ) cr i t = (Xp/ko) 1 ^ 2 , so 
that the system is always stable if c s > (c s ) cr it and becomes unstable to large wavelengths 
if c s < (c s ) cr it. Therefore, the system experiences a phase transition from a homogeneous 
distribution to an inhomogeneous distribution when the velocity of sound passes below a critical 
value (c s ) cr it. In the usual Jeans problem, the system is always unstable to large wavelengths 
since (c s ) crit = +00. The condition of instability corresponds to k < kj where kj = (Ap/c 2 ) 1 / 2 
is the Jeans wavenumber. For ko 7^ and c s < (c s ) cr i t , the condition of instability is k < k m = 
\Jk 2 — k\. Therefore, the effect of the shielding is to shift the instability to larger wavelengths 
with respect to the Jeans length. In order to measure the influence of the friction parameter 
£, we have introduced a wavenumber kd = £/(2c s ) and a dimensionless number F = (kd/kj) 2 . 
The square root of this number F 1 ! 2 ~ ^tp corresponds to the ratio between the dynamical 
time tjj ~ 1/y/pX and the friction time For F = 0, Eq. (j2J) reduces to the Euler equation 
describing a purely inertial evolution (to <C £~ x ) and for F — > +00, Eqs. (CQ)- (j21) lead to 
the generalized Smoluchowski equation (j4j) describing an overdamped evolution (£ -1 <C tjj). 
We have introduced a wavenumber k c (F) which separates, in the zone of stable wavenumbers 
(k > k m ), the case of purely exponential decay (k m < k < k c ) from the case of damped 
oscillations (k > k c ). We have also determined, in the unstable zone (k < k m ), the wavenumber 
fc* corresponding to the maximum growth rate. In the Newtonian model (no shielding) it is 
equal to fc* = corresponding to infinite wavelengths. In the Yukawa model, it corresponds to a 
finite wavelength given by Eq. (160]) . independent on the friction parameter £. For c s — > (c s ) cr i t , 
k* — ► and for c s — > 0, k* — > +00 corresponding to small wavelengths. We have found that 
the shielding length present in the chemotactic model solves many problems inherent to the 
Jeans analysis. Indeed, there is no Jeans swindle in the chemotactic problem and the maximum 
growth rate occurs for a, finite wavelength (when ko 7^ 0) instead of an infinite wavelength (when 
ko — 0). Therefore, the mathematical problem of linear dynamical stability is better posed in 
biology than in astrophysics since it avoids the Jeans swindle. 

The linear stability analysis performed in this paper gives the condition of instability (in Sec. 
HI we have expressed this condition of instability in terms of the density for different equations 
of state p(p) used in the literature) and describes the early development of the instability. 
When the condition of instability (JHj) is fulfilled, the perturbation grows until the system can 
no longer be described by equilibrium or near equilibrium equations. Therefore, the next step 
is to study the chemotactic collapse in the nonlinear regime to see the formation of patterns like 
clusters and filaments. A large number of studies in applied mathematics (see the extensive list 
of references given in the review [2] ) and physics [HJ dS, [161 EZ] have considered the overdamped 
limit of the model ([I])-© leading to the Keller-Segel model P|)-([S|), similar to the Smoluchowski- 
Poisson system (}9])- (jT0]) . For this parabolic model, chemotactic collapse leads to the formation 
of round clusters. The evolution of an individual cluster in the nonlinear regime can be studied 
by considering spherically symmetric solutions of the Keller-Segel model. The standard Keller- 
Segel model (1921) leads to the formation of Dirac peaks (for d > 2) [U Q3J [TB] . In the regularized 
model (I102p . the Dirac peaks are replaced by smooth aggregates [381 EZ]- These aggregates 
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interact with each other and lead to a coarsening process where the number of clusters N(t) 
decays with time as they collapse to each other. This process may share some analogies with 
the aggregation of vortices in two-dimensional turbulence [ID] . We expect that the decay of the 
number of clusters depends on the effective range of interaction mediated by the chemical A;^ 1 
(shielding length). If the shielding length is small the clusters do not "see" each other and the 
decay of N(t) should be slowed down or even stopped. In that case, we obtain a quasi stationary 
state made of clusters separated from each other by a distance of the order of the shielding 
length k^ 1 . If we take into account inertial effects, using the hyperbolic model ([I])-© instead of 
the parabolic model (jl])-©, the collection of isolated clusters is replaced by a network pattern 
with nodes (clusters) separated by chords (2TJ, [22], [23] • The filaments between two nodes have 
a length of the order of k^ 1 . Again, the number of nodes should decay with time. However, if 
the shielding length k^ 1 is small, the evolution is slowed down and we get a quasi equilibrium 
state with a filamentary structure corresponding to the initiation of a vasculature [2T| . 
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